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ABSTRACT 

One of the main goals of investigations using present and future giant extensive air shower (EAS) 
arrays is the mass composition of ultra- high energy cosmic rays (UHECRs). A new approach to the 
problem is presented, combining analysis of arrival directions with the statistical test of the paired 
EAS samples. An idea of the method is to search for possible correlations of UHECR masses with 
their separate sources, for instance, if there are two sources in different areas of the celestial sphere 
injecting different nuclei, but fluxes are comparable so that arrival directions are isotropic, the aim is 
to reveal a difference in the mass composition of CR fluxes. The method is based on a non-parametric 
statistical test - the Wilcoxon signed-rank routine - which does not depend on the populations fitting 
any parameterized distributions. Two particular algorithms are proposed: first, using measurements 
of the depth of EAS maximum position in the atmosphere; and second, relying on the age variance 
of air showers initiated by different primary particles. The formulated method is applied to the 
Yakutsk array data, in order to demonstrate the possibility of searching for a difference in average 
mass composition of the two UHECR sets, arriving particularly from the supergalactic plane and a 
complementary region. 

Subject headings: cosmic rays - instrumentation: miscellaneous - methods: data analysis 



1. INTRODUCTION 

The origin of ultra-high energy cosmic rays (UHECRs) 
is a long-standing challenge for astrophysics. The energy 
spectrum of particles constituting cosmic rays (CRs) is 
measured up to and slightly above 10 20 eV (= 100 EeV) 
(lAbbasi et alj|2008t ISaJamida et afll201lHEgorova et al.l 
|2004|) . but no evidence has been revealed till now, ei- 
ther of the sources or the sort(s) of highest energy par- 
ticles, owing mainly to the arrival direct i on distribution , 
which is nearly isotropic (jCroninl [20051 : iGriederl 120101 ) . 
However, some hints have been found recently of the 
possible correlation of UHECR arrival directions with 
nearb y active galactic nuclei ( AGN) at energies E > 56 
EeV (|PAO collaboration|[200l . 

Disputable estimates of the mass composition of 
the highest energy CRs were given, based on the 
average depth of EAS maximum, x m , and variance, 
a(x m ), measured by the Pierre Auger Observatory 
(PAO), the High Resolution Fly's Eye (HiRes), 
the Telescope Array (TA), and the Yakutsk array 



(lUHECR composition Working Group report at CERN Confej ] egL qg ^ n nuc i e j to 



sources of UHECRs. 

The paper is structured as follows: the new method is 
described in Section 2, where two algorithms are formu- 
lated using the depth of EAS maximum and the shower 
age measurements, and a non-parametric statistical test 
is applied to distinguish a pair of samples. In Section 3, 
the method is tested and applied to the Yakutsk array 
data observed in the energy range above 1 EeV. Conclu- 
sions are given in the final section. 

2. SEARCHING FOR A DIFFERENCE IN MASS 
COMPOSITION OF UHECRS ARRIVING FROM 
COMPLEMENTARY CELESTIAL REGIONS 

The method is aimed at the possible differences in the 
mass of UHECRs arriving from different regions in the 
celestial sphere. For example, in models of CR acceler- 
ation by shocks in AGN relativistic jets the mechanisms 
are proposed where the maximal energy of CRs is pro- 
portional to the particle charge, and wh ere protons coul d 
be accelerated up to energy ~ 100 EeV (|Berezhko l2009ft . 



l2012f) . One possible explanation of the diverging results 
may be the different average masses of UHECRs, 
observed in different fields of view of the instruments. 

In this paper, another approach is used to search for 
the possible correlation of UHECR masses with their 
sources, based on extensive air shower (EAS) observ- 
ables, namely x m and the shower age, r, varying with 
the mass of primary particles. Here, 'age' means the 
stage of the cascade development at the detector level, 
xq : t = xq sec 9 /x m , where 9 is th e inclination angle of 
the shower axis (jlvanov et al.ll20l"lf ). This method is con- 
venient for revealing different primary particles initiating 
EAS, rather than to search for excess flux from candidate 
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energy - 3000 EeV (in Cen A, 



Hondal 



(2009)). Observable particle energies should be reduced 
due to fragmentation and energy loss in the intergalac- 
tic medium and in the Galactic wind. However, in the 
Cen A case, the object is nearby (< 5 Mpc), so that the 
energy is not crucially degraded. 

By selecting EASs with arrival directions in the vicin- 
ity of Cen A in contrast with the complementary area 
where, presumably, protons dominate, o ne can reveal t he 
fraction of heavy nuclei, if the model of iHondal (|2009h is 
applicable. An obstacle is deflections of particles in mag- 
netic fields - it was shown that protons of GZK energy 
(~ 40 EeV) are deflected a few degrees coming from an y 
source within 100 Mpc (jSommers fc Westerhofll l2009h . 
Deflections of iron nuclei from Cen A would be greater 
by ~ 30%, so that the energy of particles detected should 
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be well above the GZK energy. 

In the transition region between galactic and extra- 
galactic components of CRs, the method can be used 
to verify the difference in mass composition of the two 
components comparing, for example, equatorial and po- 
lar regions in galactic coordinates. 

The method can be regarded as a n extension of the 
matte r tracer model proposed by iKoers fc Tinvakovl 
(2009]), for the case where different particles are supposed 
to be generated in UHECR sources. 

Our first task is to test the null hypothesis; that is, 
that there is no difference in mass composition between 
two regions. We have to compare two distributions of 
some EAS observables sensitive to composition of the 
primaries in the given energy range, in order to decide - 
whether there is a significant difference or not. 

If the difference exceeds experimental errors, then the 
only cause should be the mass composition in UHECR 
samples 1 . So, the next task would be to evaluate the 
most probable value and confidence interval of the mass 
difference. In this paper we focus on the first task. 

2.1. Non-parametric statistical test of data samples 

The distributions of measured EAS parameters are not 
described usually by the normal distribution, and more- 
over, have a specific form in each particular case, so that 
the general approach to the statistical test of the data 
samples is preferably non-parametric. The meaning of 
the term refers here to distribution-free methods, which 
do not rely on assumptions that the data are drawn from 
a known probability distribution. Non-parametric meth- 
ods can be used, for instance, in studying samples of EAS 
data where certain assumptions cannot be made about 
the original population. An additional advantage of the 
approach is that, even in the case when the use of para- 
metric methods is justified, non-parametric methods are 
easier to use, and leave less room for improper use and 
misunderstanding. 

In this paper, a Wilcoxon test for a pair of samples is 
used in data analysis. This test is one of the widely 
known, non-parametric signifi c ance tests (open-access 
description is given bv lLowrvl ([1998)). It is useful for 
deciding whether the two samples of observations belong 
to the same original distribution: the null hypothesis, 
Hq, is that the two samples are drawn from a single pop- 
ulation. 

Hereafter, we are going to consider pairs of matched 
samples, S and T, containing an equal number of mea- 
surements, N. In this case, the test is called the 
'Wilcoxon signed-rank test' (WSRT), in contrast to the 
case with independent samples, named the 'Wilcoxon 
rank-sum test', or Mann- Whitney U test. As an example 
of paired samples S and T in the case of cosmic rays, we 
can consider two samples of x m measurements in EAS 
detected in Summer and in Winter with the same array, 
with equal energies in pairs of events. Independent sam- 
ples are those measured in the same energy interval, but 
without equality in pairs. 

The test procedure is as follows. Excluding all pairs 
in samples where observed values are equal, Si = Tj, 
we reduce N to the number of pairs with unequal mea- 

1 We do not consider the case of different CR energy spectra 
from different celestial regions. 
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Fig. 1. — Distribution of x m . Experimental data are from Figure 
3 in llSan Luis et al.ll201ll) : 1407 events in 1 < E < 1.26 EeV, 
= 713 g/cm 2 , a(x m ) = 56 g/cm 2 . 

surements. Rank the differences \Si — Tj| in ascend- 
ing series, where rank is assigned as an item number 
in a series. The WSRT statistic, W, is then a sum 
of signed ranks 2 ; \W\ < 0.5N{N + 1). Under H 
the distribution of Wo is symmetric with W = and 
a w = y/N{N+ 1)(2N + l)/6. There are tabled values 
of the ratio zo = (Wo — 0.5) /aw with which to compare 
the resultant statistic. 

The probability P crit (z > zq) is a measure of the sig- 
nificance of deviation from the expected value under H . 
In the following, P cr n — 0.01 is assumed as the critical 
value for rejecting the null hypothesis. 

2.2. Statistical power of the Wilcoxon test 

In this section, one of the composition-sensitive EAS 
observables, i.e. the depth of shower maximum, will be 
used to find out the efficiency of the WSRT in the case 
of EAS samples. Different primary nuclei initiating EAS 
result in differ ent x m of the showe r (see, for example, re- 
cent review bv lKampert k, Ungerl (|201lD ). The question 
is: What is the minimal number of measurements needed 
to reject the null hypothesis at the significance level P cr it, 
having given difference, Ax m ? In other words, if an al- 
ternative hypothesis, H±, is true, that the two samples 
are drawn from different distributions with a given dif- 
ference in average values, Ax m , what is the probability 
to rule out the null hypothesis having samples of size N? 
Inverting the problem, we have to find the sample size 
needed to reject H at a confidence level 99%. 

It is known that, in general, the ratio of Wilcoxon test's 
efficiency to Student's t-test is 3 /it, but if the distribu- 
tions are far from the Gaussian and for large sample sizes, 
the W ilcoxon test can be consi derably more efficient than 
t-test (|Van der Waerden|[l957t ). 

In Figure [T] the distribution of x m , measured by the 
PAO Collaboration, is given (|San Luis et all 1201 1| ). It 
is undoubtedly non-Gaussian 3 . So, we will use another 
variable, xq, with normal distribution, which has the 
same average value and RMS deviation in the same en- 
ergy interval 1 < E < 1.26 EeV, in order not only to 

2 Sign is +1, if Si > T and is -1, if Si < Tj. 

3 The x 2 -deviation from the observed numbers is 168.6 for the 
normal approximation in 17 intervals, and is 1426 for the lognormal 
approximation, while less than 32 expected at the significance level 
Pcrit for good fit. 
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hypothesis, when an alternative hypothesis is true with a given 
difference, Ax m , in paired EAS samples. 



2. — Comparison of the efficiency of Wilcoxon and Student 
min is the minimal sample size required to reject the null 



determine the efficiency of the Wilcoxon test, but to com- 
pare it with a Student's t-test, known to provide an exact 
test for the equality of the means of two normal popula- 
tions with equal variances. 

Let the samples S and T of the size N be from normal 
population, 5Tg = 713 g/cm 2 , 'xf = T§ + Ax m , a(xs) = 
cj{xt) = 56 g/cm 2 . Samples are paired, so that Si and 
Tj are for showers of the fixed difference in energy of 
the primaries. Varying the corresponding difference of 
the mean depths, Ax m , one can calculate the minimal 
sample size, N m i n , needed to distinguish samples S and 
T by the two tests. 

Results of the Monte Carlo simulations are shown 
in Figure [2J The conclusion is that the ratio 
N ™?n dent l N ™n° xon is approximately 0.7 at large dif- 
ferences Ax m , and approaching 1 at minimal A.x m . It 
means that the t-test is more efficient at small sample 
sizes, i.e. N^^ ent can be reduced to ~ 70% of Wilcoxon 
test's number with equal statistical power. However, in 
our case, asymptotic efficiencies of the tests, defined as 
the limit of the efficiency as the sample size grows, are 
equal. The ratio is valid in the case of normal distribu- 
tions, where the t-test is applicable. 

2.3. A method based on the measurement of the depth 
of EAS maximum 

In the UHECR domain, the shower maximum 
is directly observable with air fluorescence detec- 
tors (FDs). Collabo rations working at HiRcs 

( Abbasi et al"][2010h. PAO dSan Luis et al.l[20Tl and TA 
( Tsunesada for the TA collaboration! 1201 ll) have mea- 
sured function of energy, to estimate the mass 
composition of CRs. In future detectors such as JEM- 
EUSO, Auger Next, etc, the fluorescence technique is 
planned to be used, including x m measurement, along 
with other developments. In all these experiments, the 
proposed method can be applied to search for correla- 
tions of UHECR masses with their sources. 

A relation between x rn and CR mass, A, is obvious 
in a superposition approximation, where an air shower, 
initiated by the nucleus of mass A, is treated as a super- 
position of A nucleon initiated showers of energy EJA . 
So, the depth of a shower maximum is (lLinslevl[l977h : 
%m = ^18 + Der ^{E/A), where Der is elongation rate, 
E is in EeV. Averaging it over the distribution of CR 
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Fig. 3. — Relation of elongation rates in EAS initiated by 
nuclei. CORSIKA simu lation results (taken from the review of 
IKampe rt & linger (20TJ)) for models are shown by three points 
calculated at E = 0.1/1/10 EeV. Superposition approximation re- 
sults in \dxm/dln E\ = |9x m /91nA|. 



energy and mass, J(E, A), we have 

x^ = xis + D E R.lnE - D ER lnA. 



(1) 

Assuming A as the slowly changing function of energy, 
we can use the relation 



= D ER A\nA 



(2) 



in a narrow energy interval. 

The accuracy of the estimations ([1]) and ^ can be 
evaluated in comparison with the model si mulation re- 
sults . Monte Carlo codes such as CORSIKA (jHeck et al.1 
1998), with implemented hadronic interaction models, 
give a more realistic description of the cascade than the 
superposition approximation. In Figure [3] a compari- 
son is presented of the results of CORSIKA simulations 
(QGS.ietll, Sibyll2 .1 and EPOSvl.99 are implemented 
|Kampert fc Ungerl |20Tl1 ) ) with equation ([J). The dif- 
ference of QGSjet results from that of equation (pj is 4% 
to 25% in the energy interval (0.1, 10) EeV; for Sibyll and 
EPOS models the divergence is not greater than 2%. 

Although the superposition approximation is not nec- 
essarily needed for our method to be applicable, it is con- 
venient to accept it for the simplicity of conclusions. So, 
we can use the results of Sibyll2.1 or EPOSvl.99 models, 
together with equations (JT]) and ([2]) within the intervals 
E e (0.1, 10) EeV, A G (1, 56), where the discrepancy is 
less or equal to 2%. 

A straightforward approach to estimation of a differ- 
ence in average mass composition is to compare x m of 
two samples of showers in the fixed energy interval. In 
the case of normal distributions with equal dispersions 
this is t-test mentioned above. If there is a statistically 
significant difference Aa; m then one can certainly con- 
clude that it is owing to the difference in average masses 
of samples, A In A. 

Actually, this approach is inefficient because of x m 
rising with energy, and mass-dependent RMS devia- 
tion a{x m ). The most stringent restriction is caused 
by energy dependence: the two samples should be 
within a narrow energy interval, where few showers 
are detected with present-day arrays. For instance, 
the maximum number of EAS events detected in the 
lg_B interval of width 0.2, where x m o bservations are 
availa ble, is 1287 for the PAO data (|San Luis et al.l 
(20119, lg£ e (18.0,18.2)), 171 for the HiRes data 
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TABLE 1 

Maximal resolution of the differences in average x m and 
UHECR mass, A In A, with Student's 4-test, using available 

DATA FROM EAS ARRAYS. 



Detector 




A;r m 


A In A 






g/cm 2 


g/cm 2 


Siby 112.1 


EPOSvl.99 


PAO 


56 


7.2 


0.29 


0.27 


HiRes 


52 


20.0 


0.81 


0.76 


TA 


52 


32.1 


1.30 


1.21 



(lAbbasi et al.l (12010ft. \eE € (18.0, 18.2)), and 6 8 for th e 
TA data (|Tsunesada for the TA collaboration! (|2011h . 
IgE e (18.6,18.8)). 

In Table 1 the estimations of the majorant resolution 
are given for models and datasets. Here, the i-test is as- 
sumed applicable and the mass dependence of x m disper- 
sion is neglected, so the real resolution should be worse. 
The maximal number of events detected \x\ A\gE = 0.2 
intervals are divided equally between samples; Ax m re- 
solvable by the i-test is calculated as in Figure [2j In 
order to estimate A In A, equation [5] is used. 

Wilcoxon test results are comparable or slightly weak. 
The experimental values of dx m / dhiE are influenced by 
the possible changes of mass composition with energy; so 
one has to use the model simulation of Der with fixed 
mass, in order to estimate the resolution. 

Assessment of the resolution is ambiguous. In this pa- 
per, we assume the resolution sufficient if 10% flux of Fe 
nuclei is resolved in the background consisting of protons. 
In this context, the maximal resolution in x rn is compa- 
rable with experimental uncertainties, the resolution of 
the average mass differences is insufficient, or hardly suf- 
ficient in the PAO case. Only future arrays with consid- 
erably larger aperture may improve the efficiency of this 
method. 

Another approach should be used, that is applicable 
to all sorts of x m distributions with energy- and mass- 
dependent parameters. WSRT is a promising method 
which can be adapted to the case. In order to use WSRT 
for the analysis of EASs arriving from different celestial 
regions, we have to compare the rank sum of x m in the 
series of events. Due to the limited number of events 
available at the highest energies of interest, we should 
extend the boundaries of energy interval from which to 
select showers. 

To do so, the paired samples can be used, in which ev- 
ery EAS from one sample has its counterpart in another 
sample with the same or closest energy. In this case, the 
two samples have the same x m if the original mass com- 
positions are identical, in spite of energy-dependent x^ 
and A. On the other hand, if there is a difference in the 
average mass of original populations, then Ax m should 
exist and is approximated by (fTJ). 

With known efficiency of WSRT in resolving Ax m 
(Figure [2]) and model simulations of Der, we can es- 
timate Aln A for the observational number of showers. 
Namely, simulation results with the CORSIKA code of 
dx m /dhxA for Sibyll (25, 24 and 23 g/cm 2 ), EPOS (27, 
26 and 26 g/cm 2 ) models at energies (1, 10 and 55 EeV) 
(|Kampert fc Ungerl |20H are used to relate Ax m and 
A In A in equation [2J 

In Table 2, the results for the PAO FD real data 



TABLE 2 

Minimal difference in average mass of UHECRs, Aln A, 
resolvable by WSRT based on the data observed by PAO 
fluorescence detectors and JEM-EUSO planned statistics. 







Aln A 




Experiment 


N oba 


Siby 112.1 


EPOSvl.99 


PAO (E > 1 EeV) 


6744 


0.13 


0.12 


PAO (E > 10 EeV) 


339 


0.61 


0.56 


JEM-EUSO (E > 55 EeV) 








nadir mode 


550 


1.07 


0.98 


tilted mode 


1800 


0.59 


0.54 



(|San Luis et al.l 1201 1[ ). and the estimated number of 
UHECRs (E > 55 EeV) to be detected using JEM- 
EUS O during 5 years o f orbiting the earth onboard the 
ISS (|Adams et al.ll201^ are shown. Duty cycle 0.19 and 
cloud impact 0.7 factors are accepted for JEM-EUSO 
exposure. The resolution in x m is assumed to be 120 
g/cm 2 . 

The number of events is only sufficient in the PAO FD 
case for the threshold energy 1 EeV. The data of JEM- 
EUSO will not be convenient for this kind of analysis, 
owing mainly to poor resolution in the depth of shower 
maximum. 

To increase the number of events under analysis at the 
highest energies, the data of the surface detectors (SD) 
can be used. This possibility will be discussed in the next 
section. 

2.4. A method based on the age variance of air 
showers initiated by different primary particles 

There are some shower parameters measured by the 
surface detectors of the EAS arrays, sensitive to the nu- 
clear composition of the primary particle: muon content, 
shower front curvature, etc. It seems that one of the most 
appropriate parameters is the shower age. It is related 
explicitly to x m A , hence, to the primary mass, and can 
be estimated in each shower using the universal relation 
with lateral distribution (LD) parameters, measured by 
the SDs, in particular the s lope, rj, of the charged parti- 
cles LD (|Ivanov et a l. 2011). Only electromagnetic com- 
ponent detectors are needed in the surface stations in 
this case, so the method can be applied to arrays with 
no muonic and other component detectors. 

As in the previous section, having two paired samples 
of EAS ages one can apply WSRT in order to decide 
whether there is an appreciable difference in the mean 
age of the showers in samples S and T, or whether these 
samples are drawn from the same distribution (the null 
hypothesis). To do so, one has to compare the ranks in 
samples of ages. Consulting with statistical tables about 
the deviation of statistic from the expected value, one 
can accept or reject the null hypothesis at the significance 
level specified. 

Using the age variance instead of Ax m in EASs initi- 
ated by different primary particles, one has an additional 
variable to fix: the zenith angle. So, the paired samples 
should be selected with close or equal energies and zenith 
angles. This can be done as follows: initially, one has two 
subsets of EASs arrived from different sources, and pre- 

4 for the fixed zenith angle 
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TABLE 3 

Estimation of the minimal difference in the shower age, 
At, and average mass, Aln A, of UHECRs observed by PAO 
surface detectors, resolvable by WSRT. 



1.0<sec(8)<l.l 
1.0<E<1.78 EeV 









Aln A 




E thr , EeV 


Nobs 


At 


Sibyll2.1 


EPOSvl.99 


3 


63376 


0.0015 


0.04 


0.03 


10 


4790 


0.0061 


0.16 


0.15 


25 


608 


0.0102 


0.28 


0.26 



sumably, with different masses. For each shower from 
the minor subset, a counterpart shower should be found 
in another subset with equal or closest 9 and E, then 
ages of these EASs should be collected in corresponding 
samples S and T. Selected showers should be removed 
from subsets. Repeating operations till the end of the 
minor subset, one assembles a pair of samples of size N 
congruous, on average, at N 3> 1, to the circumstances 
of zenith angle and energy. 

A relation between shower age and LD slope can be 
measured experimentally by PAO and TA fluorescence 
and surface detectors working together in the same show- 
ers. Meanwhile, an estimation can be used, derived by 
llvanov et aD (|2011h using the CORSIKA code with im- 
plemented SIB YLL2 . 1 /UrQMD models. 

A minimal mass difference, resolvable by WSRT ap- 
plied to age samples, is limited by uncertainty in age es- 
timation and zenith angle. Neglecting the uncertainty of 
9 measurement in PAO data (dcos9/ cos6> < 0.01), and 
using the number of EAS events detecte d during the pe- 
riod between 1.01.2004 and 31.12.2010 (jSalamida et all 
1201 If ), we have the estimation of the minimal difference 
in average mass, resolvable by WSRT applied to PAO SD 
data (Table 3). Here, the minimal difference in shower 
age resolvable by the method is used, which in turn, is 
related to the difference in x m : \At/t\ — \Ax m /x m \ for 
a fixed zenith angle. Ax m is calculated applying the test 
for a given number of EAS events in samples ~ N a i, s /2; 
average values of x m are given by t he approximation of 
the PAO data (|San Luis et alll2011l ). 

Comparison of Tables 2 and 3 shows that, in general, 
the number of events detected with PAO SDs is suffi- 
cient (contrary to FDs with reduced duty cycle) to apply 
WSRT to the shower age and mass variance of EAS pri- 
maries, in the energy range above 10 EeV, and may be 
applicable at energies above 25 EeV. Concerning future 
arrays, large apertures, sufficient resolution in x m and/or 
t, and additionally, the presence of the surface detectors 
would be essential conditions of the applicability of the 
method. 

3. APPLICATION OF THE METHOD TO 
ANALYSIS OF THE YAKUTSK ARRAY DATA 

In this section, the trial run of the method is realized 
with the Yakutsk array data. 

3.1. The Yakutsk array 

The Yakutsk array detects EAS of cosmic rays in the 
energy interval from 1 PeV to 100 EeV. The array is 
located at 61.7°N,129.4°E, 105 m above sea level (1020 
g cm -2 ). It consists of 71 surface and 6 underground 
scintillation detectors of charged particles (electrons and 




Fig. 4. — A distribution of the LD slope, r). 3058 EAS events 
detected with the Yakutsk array (points) are selected in the energy 
and zenith angle intervals, defined in the upper left corner. Normal 
approximation with the same mean value and dispersion is plotted. 

muons), and 49 detectors of the air Cherenkov light. The 
total area covered by detectors with 500 m separation is 
~ 10 km 2 . The array has been operating since 1970, and 
approximately 10 6 showers of the primary energy above 
about 10 15 eV have been detected. The highest energy 
event (E - 100 EeV) detected was 7.05.1989 with an axis 
within the array area, at zenith angle 9 = 59°. More ex- 
tended descr iption of the a r ray an d results obtai n ed ca n 
be foun d in lEgorova et all (|2004T ): llvanov et all (|2009h : 
llvanovl (12010H 

3.2. Slope of the lateral distribution function of 
charged particles in EAS 

An air shower cascading higher in the atmosphere ('old' 
shower, r > 1) has a broad and flat lateral distribution 
of secondary particles at the observational level, while a 
'young' one (r < 1) has a steep LD. It was shown pre- 
viously that the LD parameters (slope of the Cherenkov 
light and charged particles lateral distrib utions, etc.) can 
be used as indicators of the shower age dDvakonov et all 
ll^lScnmidt et alJl20M llvanov et aT1l2011ll 

In this work, the LD slope of charged particles detected 
with scintillators of the Yakutsk array is used to estimate 
the shower age in the given energy and zenith angle inter- 
vals. Cherenkov light data are not used because of the 
small sample size of showers having this kind of signal 
detected. The same reason concerns the muon detectors 
data of the array. 

The dataset used to analyze the slope parameter con- 
sists of EAS events collected during the period 1974 - 
2004, with energies from 1 to 100 EeV, zenith angles 
9 < 50° and axes within the array area. Inclined events 
beyond 50° are rejected because of substantial fraction 
of muons in the distribution of charged particles mea- 
sured. In order to estimate LD slope of each shower in a 
set, additional selection criteria were applied: i) at least 
4 stations in the core distance interval r £ (200, 1000) 
m should have particle density above a threshold; and ii) 
the slope calculated using the least square method should 
be in the interval rj € (—8,0). A total number of events 
survived after rejections is 19600. 

A distribution of slopes in the narrow interval of energy 
and zenith angle is illustrated in Figure 2] Normal ap- 
proximation is rejected by the Pearson's \ 2 test because 
of the test-statistic equal to \ 2 = 1697.8, while 23.2 is 
expected at the significance level P cr it with 10 degrees of 
freedom. 
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Fig. 5. — Efficiency of WSRT applied to paired LD slope sam- 
ples. The minimal sample size, N m i n , sufficient to reject the null 
hypothesis at P cr it when a difference in the mean slope of samples, 
Ar/, is given. 
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Fig. 6. — The probability of two LD slope samples to be drawn 
from the same distribution. 



The same procedure as in Section 2.4 is used to ap- 
ply the Wilcoxon test to paired samples, except for the 
shower age replaced by the LD slope here. Sample S 
consists of showers with arrival directions in supergalac- 
tic (S G) 'pancake' or 'dum bbell' structure around the SG 
plane (|Lahav et al.ll2000l) . Sample T consists of all other 
showers. A hypothesis tested is that UHECR sources 
in SG pancake emit nuclei, whereas from other (distant) 
sources protons arrive. Our task is to ascertain if there 
is an appreciable difference in average mass of UHECRs 
in two samples. 

In order to reveal the reliability bounds of WSRT in 
our particular case, we have applied a procedure to pair 
of EAS event samples selected in the narrow primary en- 
ergy bin E ~ 1 EeV 5 , arrival directions within the whole 
sky, but in adjacent zenith angle intervals, so that paired 
showers are of the same energy, and with given differ- 
ence in sec 6. This results in the same x m but different 
ages of paired EASs. Due to the universality of EAS, 
the average shower age is connected with the LD slope 
(jlvanov et alJl201lD . 

Shifting adjacent zenith angle intervals, we can assign 
the difference in average LD slopes of samples. Then, 
we have to determine the minimal subsample size, iV m j ra , 
sufficient to reject the null hypothesis that there is no 
difference in our source samples of slopes, by applying 
Wilcoxon test to paired subsamples. It is the efficiency 
of WSRT in the case of paired LD samples. In Figure [5] 
this limit is shown as a function of the given difference in 
average slopes. Experimental errors in charged particle 
densities measured by scintillators result in uncertainties 
of differences, Ar/, and in dispersion of points in the plot. 
The distribution of signed rank sum in the case of null 
hypothesis is calculated extracting two paired subsam- 
ples of size N rn i„ from a single sample of LD slopes. 

Although in Optical Redshift Survey and IRAS 1.2- 
Jy redshift survey data the density contrast in the local 
galaxy density field is aligned along SG axe s with radius 
40 h' 1 Mpc and thickness of 20 h' 1 Mpc (jLahav et al.l 
2000), we have used three variants of the pancake an- 
gular boundaries in supergalactic latitudes: b$G < 
15°/30°/45° in order not to miss the possible correla- 
tion of UHECRs with the SG plane proposed by IStanevi 



(120081) . 

The results are shown in Figure [51 In all cases we 
cannot reject the null hypothesis: there is no appreciable 
difference in LD slopes of EAS samples. Only at E ~ 20 
EeV, b S G < 30° is there the minimum of a probability, 
but it is greater than P cr it- 

In Table 4 the minimal differences in EAS age-related 
parameters are given, resolvable by WSRT, with the sam- 
ple sizes determined by the threshold energy. Simulation 
results concerning relations among the LD slope and the 
shower age and x „, derived with the Sibyll model, are 
used in this case (jlvanov et al.l [201l[) . A conclusion to 
be drawn is that the Yakutsk array data are too scanty 
to distinguish possible difference in mass composition of 
CRs, either at the threshold energy 10 EeV or 26 EeV, 
corresponding to C and Fe, respectively 6 . Instead, the 
data are sufficient, and can be used to search for varia- 
tions in the composition of galactic CR at energies below 
1 EeV. 

4. CONCLUSIONS 

A new method is developed, combining the analysis 
of UHECR arrival directions with the statistical test of 
paired EAS samples. It is applicable in cases when sky 
regions can be separated, where CR mass compositions 
are presumably different. The most straightforward ap- 
plication of the method lies at high energies (E > 10 
EeV) where magnetic deflections are smallest and are 
not expected to completely isotorpize the CR sky distri- 
bution. 

TABLE 4 

Estimation of the minimal difference in EAS parameters 

DETECTED WITH THE YAKUTSK ARRAY, RESOLVABLE BY WSRT. TV 
IS THE NUMBER OF CRS ARRIVED AT |&sg! < 30°; DIFFERENCES: 
IN LD SLOPE, Ar], IN THE SHOWER AGE, At, IN AVERAGE CR 

mass, A In A. 



E thr , EeV 


N 


Ar? 


Ar 


A In A 


1.0 


9180 


0.014 


0.016 


0.33 


10.0 


114 


0.120 


0.133 


3.10 


26.0 


17 


0.340 


0.378 


9.16 



5 where CRs are presumably of Galactic origin and homogeneous 
in composition as a result of confinement in magnetic fields 



6 Here, the rigidity 1 EeV is assumed as a very minimal li mit to 
charged particles not confined in the Galaxy {Horandcl 2 008T) . 
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The method is based on a non-parametric statistical 
test, WSRT, which does not depend on the populations 
fitting any parameterized distributions. Two algorithms 
are proposed: first, using measurements of the depth of 
EAS maximum position in the atmosphere; and second, 
relying on the age variance of air showers initiated by 
different primary particles. 

The efficiency of WSRT is estimated with the distribu- 
tions of shower maximum and the slope of lateral distri- 
bution in EAS. It is shown that the efficiency of Wilcoxon 
and Student tests are approximately equal in the case of 
normal distribution of the variable. 

It is also shown that the data amount, concerning x m 
measurements with present-day EAS arrays, and even 
planned for JEM-EUSO telescope, is not sufficient to 
distinguish a 10% flux of iron nuclei from the proton 
background at the significance level 0.01. However, mea- 
surements of EAS parameters related to the shower age 
with the surface detectors, specifically in the PAO exper- 
iment, provide the bulk of the data able to reveal such a 
flux. 

A trial run of the test is performed with the Yakutsk 
array data in order to search for a possible difference in 
the mass of CRs arriving from the supergalactic plane, 
and a complementary region. No significant difference is 



found in the LD slope parameter of the two subsets of 
EAS events detected with different threshold energies. 

Using the correlation between the LD slope and the 
shower age, and maximum position provided by the 
Sibyll model simulations of EAS initiated by different 
primary particles, the estimations are given of minimal 
differences in these parameters, as well as the average 
mass of CRs, resolvable by the WSRT. The number of 
EAS events detected with the Yakutsk array is found 
to be insufficient to indicate any difference in masses of 
CRs, in the energy range above 1 EeV. Instead, it can be 
used to search for variations of galactic CR composition 
below this energy. 

The method developed is quite general and can be ap- 
plied at other energies and to the various EAS parame- 
ters measured at existing and future observatories. 
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12193. 



REFERENCES 



Abbasi, R. U., et al. 2008, PRL, 100, 101101 

Salamida, F. for the PAO collaboration 2011, in Proc. 32nd Int. 

Cosmic Ray Conf. (Beijing, China: IUPAP), 2, 145 
Egorova, V.P. et al. 2004, Nucl. Phys. B (Proc. Suppl.), 136, 3 
Cronin, J. W. 2005, Nucl. Phys. B (Proc. Suppl.), 138, 465 
Grieder, P.K.F. 2010, Extensive Air Showers. High Energy 

Phenomena and Astrophysical Aspects. (2nd ed.; Berlin, 

Heidelberg: Springer- Verlag) 
PAO Collaboration 2007, Science, 318, 938 

Bcllido, J., et al., 2012, in Proc. Int. Symp. Future Directions in 

UHECR Phys. (Geneva, CERN) 7 
Ivanov, A. A., Pravdin, M.I., & Sabourov, A.V. 2011a, Int. J. Mod. 

Phys. D, 20, 1539; Ivanov, A.A., Pravdin, M.I., & Sabourov, A.V. 

2011b, in Proc. Int. Cosmic Ray Conf. (Beijing, China: IUPAP), 

2 35 

Berezhko, E. G. 2009, APJ, 698, L138 
Honda, M. 2009, APJ, 706, 1517 

Sommers, P., & Westerhoff, S. 2009, NJP, 11, 055004 

Koers, H.B.J., & Tinyakov. P. 2009. JCAP. 0904. 003 

Lowry, R,. 1998, http://vassarstats.net/textbook/ 

Kampert, K.-H., & Unger, M. 20l2, ApPh, 35, 660 

Van Der Waerden, B.L. 1957, Mathematische Statistik, (Berlin, 

Gottingen, Heidelberg: Springer- Verlag) 
San Luis, P. F. for the PAO collaboration 2011, in Proc. Int. Cosmic 

Ray Conf. (Beijing, China: IUPAP), 2, 105; PAO Collaboration 

2010, PRL, 104, 091101 



Abbasi, R. U., et al. 2010, PRL, 104, 161101; HiRes Collaboration 

2005, APJ, 622, 910 
Tsuncsada, Y., et al. 2011, in Proc. Int. Cosmic Ray Conf. (Beijing, 

China: IUPAP), 2, 246; Tameda, Y., et al. 2012, in Proc. Int. 

Symp. Future Directions in UHECR Phys. (Geneva, CERN) 75 
Linsley, J. 1977, in Proc. 15th Int. Cosmic Ray Conf. (Plovdiv, 

Bulgaria: IUPAP), 12, 56 
Heck, D., et al. 1998, Report FZKA (Forschungszentrum 

Karlsruhe), 6019 
Adams, J.H., et al. 2012, arXiv:1203.3451 

Ivanov, A.A., Knurenko, S.P., & Sleptsov, IYe. 2009, NJP, 11, 
065008 

Ivanov, A. A. 2010, APJ, 712, 746 

Dyakonov, M.N., et al. 1979, in Proc. 16th Int. Cosmic Ray Conf. 

(Kyoto, Japan: IUPAP), 8, 174 
Schmidt, F., Ave, M., Cazon, L., & Chou, A. 2008, ApPh, 29, 355 
Lahav, O., et al. 2000, MNRAS, 312, 166 
Stanev, T. 2008, arXiv:0805.1746 
Horandel, J.R. 2008, Rev. Mod. Astron. 20, 198 



